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We present a detailed study of the effects of mesh refinement boundaries on the convergence and 
stability of simulations of black hole spacetimes. We find no technical problems. In our applications 
of this technique to the evolution of puncture initial data, we demonstrate that it is possible to 
simulaneously maintain second order convergence near the puncture and extend the outer boundary 
beyond 100M, thereby approaching the asymptotically flat region in which boundary condition 
problems axe less difficult. 


I. INTRODUCTION 

Numerical relativity, which comprises the solution of 
Einstein’s equations on a computer, is an essential tool 
for understanding the behavior of strongly nonlinear dy- 
namical gravitational fields. Current grid-based formula- 
tions of numerical relativity feature ~ 17 or more coupled 
nonlinear partial differential equations that axe solved us- 
ing finite differences in 3 spatial dimensions (3-D) plus 
time. The physical systems described by these equations 
generally have a wide range of length and time scales, 
and realistic simulations are expected to require the use 
of some type of adaptive gridding in the spacetime do- 
main. 

A primary example of the type of physical system to 
be studied using numerical relativity is the final merger 
of two comparable mass inspiraling black holes, which 
is expected to be a strong source of gravitational ra- 
diation for ground-based detectors such as LIGO and 
VIRGO, as well as the space-based LISA. The individ- 
ual black hole masses Mi and M 2 set the scales for the 
binary in the source interaction region, and we can ex- 
pect both spatial and temporal changes on these scales 
as the system evolves. The binary must be evolved for 
a time t ~ 1000M, M ~ Mi -I- M 2 , starting from a sep- 
aration ~ 10M to simulate its final few orbits followed 
by the plunge and ringdown. This orbital region is sur- 
rounded by the wave zone of scale ~ 100M, where the 
outgoing signals take on a wave-like character and can be 
measured. Accomplishing realistic simulations of binary 
black hole mergers on even the most powerful computers 
clearly requires the use of variable mesh sizes over the 
spatial grid. 

Adaptive mesh refinement (AMR) was first applied 
in numerical relativity to study critical phenomena in 
scalar field collapse in 1-D [? ]; several other related 
studies have also used AMR, most recently in 2-D [? ]. 
AMR has also been used in 2-D to study the evolution 


of inhomogeneous cosmologies [? ? ]. In the area of 
black hole evolution, AMR was first applied to a simu- 
lation of a Schwarzschild black hole [? ]. Fixed mesh 
refinement (FMR) was used to evolve a short part of a 
(nonequal mass) binary black hole merger [? ] an excised 
Schwarzschild black hole in an evolving gauge [? ], and 
orbiting, equal mass black holes in a co-rotating gauge [? 
]. AMR has also been used to set binary black hole ini- 
tial data [? ? ]. The propagation of gravitational waves 
through spacetime has been carried out using AMR, first 
using a single 3-D model equation describing perturba- 
tions of a Schwarzschild black hole [? ] and later in the 
3-D Einstein equations [? ]. More recently, gravitational 
waves have been propagated across fixed mesh refinement 
boundaries, w r ith a focus on the interpolation conditions 
needed at the mesh boundaries to inhibit spurious re- 
flected waves [? ]. 

Realistic simulations of the final merger of binary black 
holes are likely to require a hierarchy of grids, using both 
FMR and AMR. The source region would have the finest 
grids, and would be surrounded by successively coarser 
grids, encompassing the orbital region and extending into 
the wave zone out to a distance > 100M. Evolving dy- 
namical gravitational fields using such a mesh refinement 
hierarchy poses a number of technical challenges. For ex- 
ample, the gravitational waves produced by the sources 
will originate as signals in the near zone and need to cross 
fixed mesh refinement boundaries to reach the wave zone. 
In addition, “Coulombic” signals that may vary with time 
but are not wavelike in character, such as are produced 
by the gravitational potential around black holes, can 
stretch across mesh boundaries. Inappropriate interpo- 
lation conditions at refinement boundaries can lead to 
spurious reflection of signals at these interfaces; c./. [? ]. 
Additional complications can arise when the grid refine- 
ment is adaptive. 

In this paper, we use the evolution of a single 
Schwarzschild black hole with FMR as a numerical lab- 
oratory. The black hole is represented as a puncture 


without excision. Using a hierarchy of fixed mesh refine- 
ments, we are able to resolve the strong field region near 
the puncture (and demonstrate the convergence of the 
solution in this region) while locating the outer bound- 
ary at > 100M. In Sec. II we describe our methodology, 
including the numerical implementation. The treatment 
of mesh refinement boundaries is discussed in Sec. III. 
Black hole evolutions with FMR are presented in Sec. IV; 
examples are given of evolutions using geodesic slicing, 
and 1 4* log slicing with zero shift. We conclude with a 
discussion in Sec. V. 


II. METHODOLOGY 
A. Basic Equations 

We use the BSSN form of the ADM equations [? ? ]. 
These equations evolve the quantities 
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written here in terms of the physical, spatial 3-metric jij 
and extrinsic curvature Kij [? ], where all indices range 


from 1 to 

3. These quantities evolve according to 
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where here and henceforth the indices of conformal quan- 
tities are raised with the conformal metric. The lapse a 
and shift specify the gauge, the full derivative notation 
d/dt = djdt — C( 3 is a partial with repect to time minus 
a Lie derivative, and the notation “TF” indicates the 
trace- free part of the expression in parenthesis. These 
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quantities are analytically subject to the physical con- 
straints 

H = R - K ab K ab + K 2 (3a) 

pi = ( ymyn_ 7 ii 7 mn )v .^ mn (3b) 

for which we recompute the physical quantities on 
the right hand sides from the evolved quantites using 
Eqs. (1). 

Note that the covariant derivatives of the lapse in 
Eq. (2b) and in Eq. (2d) are with respect to the phys- 
ical metric,, and are used here for compactness. In the 
code, this is computed 

V m V„a = d m d n a - 4:d( m 4>d n) a 

- t^ n d k a + 2y mn i kl d k (j)d l a (4) 

only in terms of conformal BSSN quantities, and the in- 
dex of the covariant derivative is raised in Eq. (2b) with 
the physical metric. The Ricci tensor in Eq. (2d) is also 
with respect to the physical metric. We compute it ac- 
cording to the decomposition 

Rij = Rij + R% (5) 

with 
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and 

R& = -2V i V j <fi-2%V k V k <t, 

+ 4Vi<fiVj<f> - 4jijV k 0V k (f> (7) 

giving the conformal and remaining pieces of the physical 
Ricci tensor. 

There are many rules of thumb in the community re- 
garding how to incorporate the constraints into the evo- 
lution equations, and in particular, when to use the in- 
dependently evolved P as opposed to recomputing the 
equivalent quantity from the evolved metric. We have 
made our choices manifest in the writing of the equations 
here; we largely follow the rules set out in [? ]. 

B. Numerical Implementation 

For the spatial discretization of Eqs, (2), we take the 
data to be defined at the centers of the spatial grid cells 
and use standard O(Ax) 2 centered spatial differences [? 
]. To advance this system in time, we use the iterated 
Crank-Nicholson method with 2 iterations [? ], which 
gives O(At) 2 accuracy. For the evolutions with geodesic 
slicing, an analytic solution is available; this is used to 
set the boundary conditions at the outer edge of the 
grid for these runs. Evolutions using the 1 + log slic- 
ing condition employ interpolated Sommerfeld outgoing 



3 


wave conditions at the outer boundary. Overall, the code 
is second-order convergent; specific examples of this are 
given in Sec. IV below. 

We explicitly enforce the algebraic constraints that 
is trace-free and that 7=1. We enforce the trace-free 
condition by replacing the evolved variable with 

Aij — > Aij — — j mn A rnn jij 

and we enforce the unit determinant condition by replac- 
ing the evolved metric with 

i ij ->■ 7~ 1/3 7 ij 

after each ICN interation. Both of these constraints axe 
enforced in all of the runs presented below. 

We use the Paramesh package [? ] to implement 
both mesh refinement and parallelization in our code. 
Paramesh works on logically Cartesian, or structured, 
grids and carries out the mesh refinement on grid blocks. 
When refinement is needed, the grid blocks needing re- 
finement are bisected in each coordinate direction, similar 
to the technique of Ref. [? ]. 

All grid blocks have the same logical structure, with 
n x zones in the z-direction, and similarly for n y and 
n z . Thus, refinement of a 1-D block produces two child 
blocks, each having n x n y n z zones but with zone sizes a 
factor of two smaller than in the parent block. Refine- 
ment can continue as needed on the child blocks, with the 
restriction that the grid spacing can change only by a fac- 
tor of two, or one refinement level, at any location in the 
spatial domain. Each grid block is surrounded by a num- 
ber of guard cell layers that are used to calculate finite 
difference spatial derivatives near the block’s boundary. 
These guard cells are filled using data from the interior 
cells of the given block and the adjacent block. 

Paramesh can be used in applications requiring FMR, 
AMR, or a combination of these. The package takes care 
of creating the grid blocks, as well as building and main- 
taining the data structures needed to track the spatial 
relationships between blocks. Paramesh handles all inter- 
block communications and keeps track of physical bound- 
aries on which particular conditions are set, guarantee- 
ing that the child blocks inherit this information from 
the parent blocks. In a parallel environment, Paramesh 
distributes the blocks among the available processors to 
achieve load balance, minimize inter-processor communi- 
cations, and maximize block locality. 

This scheme, it is worth noting, provides excellent com- 
putational scalabality. Equipped with Paramesh, the 
scalabality of our code has been tested for up to 256 pro- 
cessors and has demonstrated a consistently good scaling 
factor for both the unigrid and the FMR runs. For uni- 
grid runs, we started with a uniform Cartesian grid of a 
certain number of grid cells, a fixed number of timesteps, 
and a certain number of PEs (Processing Elements), and 
then increased the number of PEs to run a larger job 
while the number of grid cells per PE remained con- 
stant. In this situation, we expect that the total CPU 


time taken to run the code should have a scaling factor 
of 1 under perfect conditions. In reality, communica- 
tion overhead makes the scalability less than perfect. We 
define the scaling factor to be the time expected with 
perfect scaling divided by the actual time taken. Using 
FMR, we ran the same simulations as in the unigrid case 
except that a quarter of the computational domain was 
covered by a mesh with twice the resolution. Despite the 
more complicated communication patterns, scalability in 
the FMR runs is comparable to the unigrid runs. The 
scaling factor of our code is ~ .92 for unigrid runs and 
~ .90 for FMR runs. 

For the work described in this paper, we are using 
FMR. For simplicity, we use the same timestep, chosen 
for stability on the finest grid, over the entire computa- 
tional domain; c.f. [? ]. At the mesh refinement bound- 
aries, we use a single layer of guard cells. Special atten- 
tion is paid to the restriction (transfer of data from fine to 
coarse grids) and prolongation (coarse to fine) operations 
used to set the data in these guard cells, as discussed in 
the next section. 


III. TREATMENT OF REFINEMENT 
BOUNDARIES 

A. Guard Cell Filling 

As a general rule, guard cells at mesh refinement 
boundaries should be filled to third order accuracy for 
an evolution code such as ours to maintain overall second 
order accuracy. 1 The current version of our code uses a 
third order guard cell filling scheme that is now included 
with the standard Paramesh package. This guard cell fill- 
ing proceeds in three steps. The first step is a restriction 
operation in which interior fine grid cells are used to fill 
the coarse grid guard cells. This is depicted for the case 
of two spatial dimensions in the left side of Fig. 1. The 
restriction proceeds as a succession of one-dimensional 
quadratic interpolations, and is accurate to third order 
in the grid spacing. Note that the 3-cell-wide fine grid 
stencil used for this step (nine darker points in the fig- 
ure) cannot be centered on the guardcell (grey square). 
In each dimension the stencil includes two fine grid cells 
on one side of the guard cell and one fine grid cell on 
the other. The stencil is always positioned so that its 
center is shifted toward the center of the block (assumed 
in the figure to be toward the upper left). This ensures 


1 In our terminology the “order of accuracy” refers to the order 
of errors in the grid spacing. Thus, third order accuracy for 
guard cell filling means that the guard cell values have errors of 
order Ax 3 , where Ax is the (fine) grid spacing. (Note that third 
order accurate guard cell filling was termed “quadratic” guard 
cell filling in Ref. [? ].) Second order accuracy for the evolution 
code means that, after a finite evolution time, the field variables 
have errors of order Ax 2 . 
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that only interior fine grid points, and no fine grid guard 
cells, are used in this first guard cell filling step. For the 
second step, the fine grid guard cells are filled by prolon- 
gation from the coarse grid using the stencil shown in the 
right side of Fig. 1. The prolongation operation proceeds 
as a succession of one-dimensional quadratic interpola- 
tions, and is third order accurate. In this case, the coarse 
grid stencil includes a layer of guard cells (grey squares), 
which have just been filled in the first step. 

After these first two steps the coarse and fine grid 
guard cells are filled to third order accuracy. However, 
order of accuracy is not our only concern. In practice, 
we find that this two-step guard cell filling scheme leads 
to unacceptably large reflection and transmission errors 
for waves passing through the interface. These errors can 
be suppressed by adding a third step to guard cell filling, 
“derivative matching” at the interface. With derivative 
matching the coarse grid guard cell values are recom- 
puted so that the first derivatives at the interface, as 
computed on the coarse grid, match the first derivatives 
at the interface as computed on the fine grid. The first 
derivative on the coarse grid is obtained from standard 
second order differencing using a guard cell and its neigh- 
bor across the interface. The first derivatives on the fine 
grid are computed using guard cells and their neighbors 
across the interface, appropriately averaged to align with 
the coarse grid cell centers. This third guard cell filling 
step maintains third order accuracy and greatly improves 
the transmission and reflection errors at the interface. 

B. Error Analysis 

Why should third order guard cell filling be adequate 
to maintain overall second order accuracy? This is a non- 
trivial question, and there are certain subtleties that arise 
in our black hole evolutions. To help pave the way for a 
better understanding of these issues, let us consider the 
simple example of a scalar field in one spatial dimension, 

ip — 7r (8a) 

7 r = ip" . (8b) 

Here, the dot denotes a time derivative and primes de- 
note space derivatives. These equations can be solved 
numerically using the twice-iterated Crank-Nicholson al- 
gorithm. The fields at timestep n + 1 are given in terms 
of the fields at timestep n by 

At 2 At 3 

V>” +1 = </>" + At tt? + (9a) 

A t 2 At 3 

7 r; +1 - 7T? + AtD 2 ^ + — D 2 tt? + — . D 4 ^ n (9b) 

where D 2 is a finite difference operator approximating 
the second spatial derivative. Consider for the moment a 
uniform spatial grid. If D 2 is the usual second order ac- 
curate centered difference operator, the dominant source 


of error for t/^ +1 comes from the term proportional to 
At 3 . This term has the “wrong” numerical coefficient 
as compared to the Taylor series expansion of the exact 
solution. The dominant sources of error for 7r? +1 come 
from the term proportional to At 3 , which also has the 
“wrong” numerical coefficient, and from the second or- 
der error in D 2 ip™. For a uniform grid the dominant 
error in D 2 ip^ is D A ip^ Ax 2 / 12, so the leading errors for 
a single timestep are 

= ^A* 3 D 2 7r; (10a) 

W +1 )err = ~ At (At 2 + A X 2 )D 4 ^ . (10b) 

For At ~ Ax, each of these one-time-step errors is pro- 
portional to Ax 3 . If we evolve the initial data to a finite 
time T, the (9(Ax 3 ) errors accumulate over N = T/ At 
timesteps resulting in second order errors. 2 Thus, the 
basic variables ip and i r are second order convergent on a 
uniform grid. 

On a non-uniform grid, guard cell filling introduces 
errors of order Ax 3 in ip at grid points adjacent to the 
boundary. This leads to errors of order Ax in D 2 ip ? 
and 1/Ax in D^ip 7 -. From Eq. (9b) we see that in one 
timestep i r can aquire errors of order Ax 2 . The concern 
is that these errors might accumulate over N — T/At 
timesteps to yield first order errors. This, in fact, does 
not happen. Simple numerical experiments show that ip 
and 7 r are second order convergent on a non-uniform grid 
with third order guard cell filling. We can understand 
this result with the following heuristic reasoning. The 
numerical algorithm of Eq. (9), as does any mathemat- 
ically sound numerical scheme, approximates the exact 
solution of the scalar field equations (Eq. (8)) in which 
the fields propagate along the light cone. The “bulk” er- 
rors displayed in Eq. (10) accumulate along the past light 
cone to produce an overall error of order N Ax 3 ~ Ax 2 at 
each spacetime point. Errors in guard cell filling, which 
occur at a fixed spatial location, do not accumulate over 
multiple timesteps since the past light cone of a given 
spacetime point will cross the interface (typically) no 
more than once. 

The characteristic fields for the system of Eq. (8) are 
7 r-t -ip 1 and 7 r — ip'. Numerically, the field ip' is obtained by 
evolving the ip, 7 r system for a given number of timesteps 
then taking the second order accurate numerical first 
derivative of ip. A more complete analysis shows that 
tp f , like 7 r, is subject to one-time-step errors of order 
Ax 2 due to guard cell filling. Since ip f takes its value 
from data along the past light cone, these errors do not 


2 This is a simplification. The dominant error after N timesteps 
includes other terms of order Ax 2 in addition to the product of 
N and the one-time-step error. These other terms include, for 
example, the product of an order TV 3 coefficient and a one-time- 
step error of order Ax 5 . 
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FIG. 1: Guard cell filling in two spatial dimensions. In these pictures, the thick vertical line represents a refinement boundary. 
The picture on the left shows the first step, in which one of the coarse grid guard cells (grey square) is filled using quadratic 
interpolation across nine interior fine grid cells (black circles). The other coarse grid guard cells are filled using corresponding 
stencils of nine interior fine grid cells. The picture on the right shows the second step in which two fine grid guard cells (grey 
circles) axe filled using quadratic interpolation across nine coarse grid values (squares). These coarse grid values include one 
layer of guard cells. The final step in guard cell filling is to use “derivative matching” to refill the coarse grid guard cells. The 
asymmetry in the left panel is drawn with the assumption that the fine block’s center is toward the top-left of the panel. 


accumulate over the course of N timesteps and the error 
in xp 1 remains second order. 3 If, after a finite evolution 
time, we compute xp” (which is the first derivative of xp’), 
then the third order guard cell filling errors in ip' will 
yield second order errors in xp” . Thus, we can evolve the 
system of Eq. (8) for finite time on a grid with mesh re- 
finement, numerically compute the second derivative of 
xp y and find that xp” is second order accurate. 

The BSSN equations, Eq. (2c) and Eq. (2d), are similar 
to the scalar field equations, Eq. (8), with 7^ playing the 
role of xp and —Aij playing the role of 7 r. This feature was 
one of the original motivations behind the BSSN system. 
Note that the term analogous to xp” in the tt equation 
is the term contained in the trace-free part of 

the Ricci tensor, which appears on the right-hand side 
of Eq. (2d). Obviously there are many other terms that 
appear on the right-hand side of the dAij jdt equation. 
We can model the affect of these terms by including a 
fixed function on the right-hand side of the tt equation: 

xp = 7r (11a) 

7T = Xp” — X * (Hb) 

We have written the fixed function as the second deriva- 
tive of x • For simplicity we choose x to depend on x 
only, the most relevant dependence for our consideration 
of behavior across spatial resolution interfaces. The gen- 


3 The one-time-step errors for xp due to guard cell filling are order 
Ax 3 . These errors do accumulate over N timesteps to yield errors 
of order N Ax 3 ~ Ax 2 . 


eral solution of this system is then 

rp(t,x) = xp(t,x) + x{x) (12a) 

7 r(t,x) = ff (t,x) (12b) 

where xp, ft is a solution of the homogeneous wave equa- 
tion (Eq. (8)). 

The extended model system of Eq. (11) can be solved 
numerically with the discretization 

^ n+1 = iPJ + At^ + ... (13a) 

7 +1 = 7 + At (D 2 ^ - D 2 Xj) + ■■■ (13b) 

The higher order terms in At, not show here, come from 
the iterations in our iterated Crank-Nicholson algorithm. 
It is important to recognize that the x" term is expressed 
as the numerical second derivative of \j stnd not as the 
discretization of the analytical second derivative, {x n )j' 
The reason for this choice is that D 2 \j mimics the effect 
of the extra terms on the right-hand side of Eq. (2d) 
which, in our BSSN code, depend on the discrete first 
and second derivatives of the BSSN variables <p and f\ 
From the discussion of the wave equation (Eq. (8)) we 
can anticipate the results of numerical experiments with 
the model system (Eq. (11)) on a non-uniform grid. For 
arbitrary initial data xp !•, the numerical solution is 

given by 

</>; = + (14a) 

7T? = 7tJ (14b) 

where xp f tj is the numerical solution of the homoge- 
neous wave equation (Eq. (8)) with initial data xp® - Xj> 
Or®. The order of convergence for xpj is determined by how 
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FIG. 2: Convergence tests for ip and ip f . The mesh interface 
is at x = 0, with the fine grid on the left and coarse grid on 
the right. The errors in ip and ip f are divided by factors of 4 
and 16 for the middle and low resolution cases. 

rapidly, as Ax -4 0, the numerical solution in Eq. (14a) 
approaches the exact solution ip(t,x) = ip(t,x) + x( x )- 
Since Xj is simply the projection of the analytic function 
x(x) onto the numerical grid, the term Xj in the numer- 
ical solution (Eq. (14a)) does not contribute any error. 
We have already determined that on a non-uniform grid 
'ip™ approaches ip(t, x) with second order accuracy. Thus, 
we expect ip ™ to be second order convergent. 

What about derivatives of ip? The order of conver- 
gence for D l ip™ is found by comparing the discrete deriva- 
tive D 1 ip™ — D 1 ^ + D l Xj to the analytic solution 
ip* = ip* + x 7 - Again, as we have discussed, D l ip™ ap- 
proaches ip' with second order errors. It is also easy to see 
that the numerical derivative D l Xj approaches x! with 
second order accuracy. Away from grid interfaces this is 
obviously true, assuming that D l is the standard second 
order accurate centered difference operator. For points 
adjacent to a grid interface, guard cell values for ip are 
filled with third order errors. These errors lead to second 
order errors in D lf ip™. Overall then, we expect second 
order convergence for D l ip ™. 

The expected convergence rates for ip and ip' are con- 
firmed by the results shown in Fig. 2. For these numerical 
tests, we chose x( x ) = exp((x - 50) / 10) and initial data 

tp(0,x) = 100e _(x+10)2/40 ° + e {x - 50)/1 ° (15a) 

tt(0,x) = I(x + 10)e- (a+lo)2/400 . (15b) 

Each set of curves shows the errors at three different reso- 
lutions, Ax = 5/16, 5/32, and 5/64, where Ax is the fine 
grid spacing. The evolution time is 20.83, corresponding 
to 200, 400, or 800 timesteps (depending on the resolu- 
tion) and a Courant factor of 1/3. 

The order of convergence for the second derivative of 
ip is determined from a comparison of D 2 ip™ = D 2 ip™ -F 
D 2 Xj and the analytic solution ip” = ip” + x”- We have 
seen that D 2 ip ™ approaches ip” with second order ac- 
curacy. The situation for D 2 Xj> however, is somewhat 
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FIG. 3: Convergence test for ip”. The solid curve is the error 
in ip” for the highest resolution run. As in Fig. 2, the interface 
is at x = 0 and the errors for the middle and low resolution 
cases are divided by factors of 4 and 16. The grid points 
adjacent to the interface do not coincide because ip” is only 
first order convergent at these points. 


different. Away from any grid interface D 2 Xj will ap- 
proach x" with second order accuracy, assuming D 2 is 
the standard second order accurate finite difference op- 
erator. But for points adjacent to the interface, and only 
those points, guard cell filling errors of order Ax 3 in x 
will lead to first order errors in D 2 Xj • Thus, we expect 
to find second order convergence for D 2 ip ™ at all points 
except those points adjacent to the interface. Points ad- 
jacent to the interface should be first order convergent. 

Fig. 3 shows the results of our convergence test for 
ip”. The spikes at the interface (x = 0) appear because 
the two grid points adjacent to the interface are only first 
order convergent. Elsewhere, the plot shows second order 
convergence. 4 

The behavior demonstrated in Fig. 3 also occurs in 
the BSSN system when we examine the convergence plot 
for the Hamiltonian constraint, Fig. 5. In graphing the 
Hamiltonian constraint HC, we are comparing a combi- 
nation of grid functions that includes second derivatives 
of the BSSN variables to the exact analytical solution of 
HC, namely, zero. We therefore expect spikes to appear 
in the convergence plot for the Hamiltonian constraint, 
and indeed they do. 

We wish to emphasize that the lack of second order 
convergence for second spatial derivatives at grid points 
adjacent to the interfaces is not due to any error in our 
code, or shortcoming of the numerical algorithm. Since 
the undifferentiated variables are second order conver- 
gent everywhere, we can always assure second order con- 


4 Although our primary concern here is in ip and its derivatives, 
it is interesting to note that 7r and 7 r' are second order conver- 
gent everywhere. The second derivative 7r" contains “packets” 
of high frequency, first order error that propagate away from the 
interface. 
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vergence of their derivatives by using suitable finite dif- 
ference stencils. For example, in computing D 2 ^ from 
xjjj we can use a second order accurate one-sided operator 
D 2 that avoids using guard cell values altogether. With 
such a choice the spikes in Fig. 3 disappear, and D 2 ^ is 
everywhere second order convergent. In our BSSN code, 
it is most convenient to compute the Hamiltonian con- 
straint using the same centered difference operator D 2 
that we use for the evolution equations. As a conse- 
quence, spikes appear at the grid interfaces in the con- 
vergence plot (Fig. 5). 

To summarize, we use a third order accurate algo- 
rithm to fill guard cells for the BSSN variables (1). As a 
consequence, second spatial derivatives of these variables 
aquire first order errors at grid points adjacent to mesh 
refinement boundaries. These first order errors show up 
as spikes in a convergence plot for quantities that depend 
on second spatial derivatives, such as the Hamiltonian 
constraint. The key result or our analysis is a demon- 
stration that the first order errors in second derivatives 
do not spoil the overall second order convergence of the 
evolved variables (1), in spite of the fact that second spa- 
tial derivatives appear on the right-hand sides of the evo- 
lution equations (2). 


IV. BLACK HOLE EVOLUTIONS 

Black hole spacetimes are a particularly challenging 
subject for numerical study. Practical applications will 
require that our FMR implementation perform robustly 
under the adverse conditions which arise in black hole 
simulations. In this section we will test that our tech- 
niques perform convergently and accurately in the pres- 
ence of strong field dynamics and singular “punctures” 
associated with black hole evolutions and that these 
methods can be stable on the timescales required for in- 
teresting simulations. 

The puncture approach to black hole spacetimes gener- 
alizes the Brill-Lindquist [? ] prescription of initial data 
for resting black holes. In this approach, the spacetime 
is sliced in such a way as to avoid intersecting the black 
hole singularity, and the spatial slices are topologically 
isomorphic to M 3 minus one point, a puncture, for each 
hole. The punctures represent an inner asymptotic re- 
gion of the slice which can be conformally transformed 
to data which are regular on JR 3 . In this way a rest- 
ing black hole of mass M located at r = 0 is expressed 
in isotropic spatial coordinates by 7 with 
factor ^bl = 1 + M/2r and Kij = 0. A direct gener- 
alization of this expression for the conformal factor can 
be used to represent multiple black hole punctures, and 
data for spinning and moving black holes can be con- 
structed according to the Bowen- York [? ] prescription. 
A key characteristic which makes this representation ap- 
propriate for spacetime simulations is that with suitable 
conditions on the regularity of the lapse and shift, the 
evolution equations imply that time derivatives of the 


data at the puncture are regular everywhere despite the 
blow up in ^bl at the puncture. Numerically, we treat 
the punctures by the prescription similar to that given 
in [? ]. In the BSSN formulation, this amounts to a 
splitting of the conformal factor exp (40) into a regular 
part, exp(40 r ), and non-evolving singular part, exp(40 s ), 
given by ^bl- The numerical grid is staggered to make 
sure the puncture does not fall directly on a grid point, 
and to avoid the large finite differencing error, the deriva- 
tives of <f) s are specified analytically. 

We study two test problems, each representing a 
Schwarzschild black hole in a different coordinate sys- 
tem. The first case, described in Sec. IV A, is a black 
hole in geodesic coordinates, in which the data evolve as 
the slice quickly advances into the singularity. In this 
case we have an exact solution for the development of 
the spacetime with which we can compare for a direct 
test of the simulation. Our next test, in Sec. IV B is in- 
tended to study the performance of our FMR interfaces 
under the condition that strong-field spacetime features 
pass through the interfaces. In this test we use a variant 
of the “1+log” slicing condition to define the lapse, a, 
with vanishing shift, /?* = 0. This gauge choice allows 
the slice to avoid running into the singularity, but causes 
the black hole to appear to grow in coordinate space so 
that the horizon region passes though our FMR interface. 


A. Geodesic Slicing Results 

1. Analytic Solution 

Geodesic coordinates are obtained by setting unit 
lapse, a = 1, and vanishing shift, /3 l — 0. This im- 
plies that the gridpoints will follow geodesic trajectories 
through the physical spacetime. We obtain the analytic 
solution in terms of these geodesics. 

The Schwarzschild geometry in standard coordinates 
is given by 

ds 2 — ~a 2 dt 2 4- a~ 2 dR 2 -j- R 2 dfl 2 , (16) 
where a 2 = (1 - 2 M/R). 

To express this metric in geodesic coordinates, consider 
a spatial Cauchy surface E 0 in a 4-manifold M and a 
congruence of radial geodesics crossing E 0 . Every point 
(Ro,Q 0 ) in S 0 can be considered the “initial” point of 
a geodesic. Let the affine parameter r for each geodesic 
be zero at E 0 . Considering subsequent slices of constant 
proper time, we can set a global time r over M and 
foliate it with the slices E r . 

Since each (Rq,£Iq) labels a unique geodesic, which 
in turn is parameterized by r, Rq can be considered a 
new radial coordinate that naturally pairs with the time 
coordinate r. To emphasize Ro's role as a free coordinate, 
we rename it p. 

Now we will derive the metric components in this coor- 
dinate system. The normalization of the tangent vector, 
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n a = ( djdr) a implies that g rr = — 1. The g rp metric 
component is the p component of the shift vector. By 
means of the geodesic equation n a V a n b = 0 we can show 
that g TPiT = 0. As initially /?\ = 0, we conclude that 

9rp — 0. 

Together with Eq. 16 and the above considerations we 
have that 


[? ? ], the extrinsic curvature can be viewed as the rate 
of change of the spatial metric 


K ab = - 


1 dy ab 

2a dr 


( 22 ) 


Using the previous expression with unit lapse and zero 
shift we obtain: 


_ jdR/dpf gRR 
y " (dt/dr) 2 g tt 


(17) 






(23) 


The final ingredient is the conserved energy associated 
to each geodesic. As the Schwartzchild geometry has a 
timelike Killing field = du we have that the energy 
E = = — {dt/dr) gtt is conserved along the 

trajectory. Initially we have E — — where = 
9ab |t=o* Therefore: 


9pp — 9rr 



To proceed, note that if we have a function / = 
/(it, u, w) = 0 defining u as an implicit function of v and 
w, we can use the chain rule and the implicit function 
theorem to show that 


du 

dv 

d 2 u 

dwdv 


df/dv 
’ df/du 

d£ 

du 


d 2 f 

dvdw 


+ 


d 2 fd f df 
du 2 dv dw 


(24) 



The trajectories R(r,p) can be explicilty computed 
from the n 2 = — 1 and the enrgy conservation for the 
above geodesics, and we obtain [? ? ]: 


oV 2 


(2M) 1 / 2 


arccos * 


1-MJ) 


0. (19) 


This expression provides an implicit definition for R — 
R(p, r), which is easily inverted numerically to high pre- 
cision. 

To perform numerical evolutions the geodesic coordi- 
nates have a drawback: the physical singularity is already 
present on the initial slice (r — 0) at p = 0. We can avoid 
this problem by going to isotropic coordinates by means 
of the transformation p = r (1 + M/2r) 2 . We see that 
both r = 0 and r — oo map to p — oo, and for real 
r the minimum value of p is p — 2 M (the horizon) at 
r = M/2, now the surface closest to the physical singu- 
larity on the initial slice. Substituting into Eq. (19) we 
see that geodesics originating on this surface reach the 
physical singularity, R = 0, at time r = 7 rM, defining 
the maximum temporal extent of our coordinate system. 

Returning to the metric, the transformation to 
isotropic coordinates gives us 



From Eq. (20) we can express the final metric as: 

ds 2 = - dr 2 + + + R 2 d£l 2 . (21) 

Expressions for the extrinsic curvature, which have not 
previously appeared in the literature, can be derived in 
a similar manner. As we know from the ADM formalism 


(df) ' 2 (&f_ d 2 f df\ 

\du) \dvdu dw dudw dv ) ' 


(25) 


Taking for / the left hand side of Eq. (19), and, noting 
that df /dr = 1, we conclude that: 


m \ 4 a/ (diy 3 

2 r ) dp V dR ) 

x d lL Pfdffdfy 1 
dpdR dR? dp \dRJ 

and 



K » - r (I) 


K" = K" sin 2 (0). 


(26) 


(27) 


There are no off-diagonal terms. Observe that we have 
only partial derivatives of / that can be obtained analyt- 
ically from Eq. (19), and easily evaluated numerically. 


2. Numerical Results 

Now we turn to the numerical evolution of a single 
puncture black hole. As explained in the previous sec- 
tion, at t = 7 rM the slice will reach the physical sin- 
gularity and because we are not performing any excision, 
this sets the maximum duration of our evolution. Never- 
theless, t w 3 M is enough to test the relevant features of 
FMR evolution and provide us with a simple analytical 
solution. 

In these simulations we locate the puncture black hole 
at the origin. The spherical symmetry of the problem 
allows us to restrict the simulation to an octant do- 
main with symmetry boundary conditions, thereby sav- 
ing memory and speeding up the evolutions. The outer 
boundaries are planes of constant x, y or z at 128M each. 
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FIG. 4: Evolution of the conformal metric component 7 XX , 
for a geodesically sliced puncture, shown at t = 0.5, 1.0, 1.5, 
2.0, and 2.5M. 


As noted before, one of the important features of FMR 
is to enable the outer boundary to be very far, of the 
order of 100M. In this case we can apply our exact solu- 
tion as the outer boundary condition, though this distant 
boundary is completely irrelevant for the most interesting 
strong-field region. 

In this test we are mainly interested in how the FMR 
boundaries behave near the puncture and under strong 
gravitational fields. Even with the outer boundary far 
away we can, by applying multiple nested refinement re- 
gions, highly resolve the region near the puncture, as 
is required to demonstrate numerical convergence. To 
achieve the desired resolution near the puncture, we use 8 
cubical refinement levels, locating the refinement bound- 
aries at 64M, 32 M, 16M, 8M, 4 M, 2 M and 1 M. To 
test convergence we will examine the results of three 
runs with identical FMR grid structures, but different 
resolutions. The lowest resolution run has gridpoints 
A Xf = Ay/ = A Zf = M/16 apart in the finest re- 
finement region near the puncture. The medium res- 
olution run has double the resolution of the first run 
in each refinement region, and the the highest resolu- 
tion run is doubled again for a maximum resolution of 
A Xf = M/64. The memory demand and computational 
load per timestep are similar to unigrid runs of 32 3 , 64 3 
and 128 3 gridpoints. 

In Fig. 4 we plot the conformal metric component % x 
along the x-axis at times t — 0.5, 1.0, 1.5, 2.0, and 2.5M. 
The metric grows sharply near 0.5M blowing up com- 
pletely as we approach t — i rM when, in these coordi- 
nates, the r = .5M surface reaches the physical singular- 
ity. In the present context though, we are not so much 
interested in the field values of this well-studied space- 
time, as in the difference between analytical and numeri- 
cal results, to directly measure the error in our numerical 
simulation. At late times, the error is, not surprisingly, 
dominated by finite differencing error in the vicinity of 
the developing singularity. 

The plots in Fig. 5 compare the errors at t — 2.5 M 


for the three varied resolution runs described above, 
demonstrating the convergence of 7 XX , A xx , the Hamil- 
tonian constraint, and the x-component of the momen- 
tum constraint. The error curve for the medium and 
low resolution runs have been divided by 4 and 16, re- 
spectively. That the curves shown lie nearly atop one 
another is an indication of second order convergence, t.e. 
that the lowest order error term depends quadratically 
on the gridspacing, Ax. That the remaining difference 
between the adjusted curves near x = 0.5M seems also 
to decrease quadratically is an indication that the next 
significant error term is of order Ax 4 . We achieve con- 
vergence to the analytic solution everywhere, from very 
near the puncture, through the peak region which is ap- 
proaching the singularity, and in the weak field region. 

Our particular interest is in the region near the refine- 
ment interfaces. The close-up shown in Fig. 7 demon- 
strates that 7 XX converges quadratically to the exact so- 
lution. Our computation of the Hamiltonian constraint, 
however, requires second derivatives and, at the point 
nearest the refinement boundary, depends on the guard- 
cells in such a way that the result is only 1st order con- 
vergent there. A close-up of the Hamiltonian constraint 
near a refinement boundary (Fig. ??) shows second or- 
der convergence at all non-adjacent points. This is as ex- 
pected according to the discussion in Sec. Ill (c/. Fig. 3). 
We have verified for these simulations that the data for 
all BSSN variables do converge at second order at the 
refinement boundaries. 


B. 1 + log slicing 

Having rigorously tested the code against analytic so- 
lutions, we now use a different coordinate condition to 
show a longer-lived run with nontrivial, nonlinear dy- 
namical behavior in the region of FMR interface bound- 
aries. For this purpose, we again use zero shift but with 
a modified “1-i-log” slicing condition given by 

/J/y 

— = -2 a<k% L K (28) 

where insertion of the factor ^bl = 1 + M/2r, originally 
recommended by [? ], improves stability as well as con- 
vergence near the puncture. With this gauge condition, 
the code runs up to t ~ 35M. For t > 35 M, the growth 
in the metric variables becomes too steep for the code to 
generate meaningful data. 

The grid set-up is the same as in the geodesic case. 
Fig. ?? shows the second order convergence of the Hamil- 
tonian constraint. The convergence of the Hamiltonian 
constraint is second order everywhere except at the re- 
finement boundaries atx~lM,x~2M where the 
convergence is only of the first order and at x ~ 0.25M 
where the code does not converge at all. The code is 
second order convergent in other refinement boundaries, 
x ~ 4M, 8M, etc. At x ~ 0.25M, the metric compo- 
nent is second order convergent (Fig. ??), but the lapse 
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FIG. 5: Convergence of the errors (numerical values minus analytic values) in ^ xx , A XX} the Hamiltonian constraint H, and 
the momentum constraint P x , for a geodesically sliced puncture, all at the time t = 2.5 M. 



FIG. 6: Hamiltonian Constraint 


(Fig. ??), exhibit high frequency noise. 

This noise requires some explanation; however, it does 
not come from the refinement boundaries. Rather it is 
generated right on the spot at x ~ 0.25 M. It becomes 
most evident at t > 10 M, first appearing in the lapse and 
K , which are directly coupled with each other, and then 
eventually in all the extrinsic curvature variables. For 
the duration of the evolutions the noise does not prop- 
agate at all. The location of the noise is always within 
x ~ 0.2M and x ~ QAM, independent of resolution and 
of the position of the FMR boundaries. The noise ap- 


pears to represent a gauge shock associated with the col- 
lapse of the lapse, and is unique to 1+log slicing with 
zero shift. Outside this region from 0.2M to QAM , all 
basic variables are second order convergent everywhere, 
including at FMR boundaries. 


V. DISCUSSION AND CONCLUSIONS 

This paper demonstrates that fixed mesh refinement 
boundaries can be located in the strong field region of a 
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FIG. 7: Close-up of the convergence of the error (numerical 
value minus analytic value) in j xx , for a geodesically sliced 
puncture, at t = 2.5 M, in the vicinity of the refinement 
boundary at x = 1M. Guardcell points are also shown. 

black hole spacetime when the interface condition is han- 
dled properly, and that using mesh refinement technol- 
ogy, therefore, is a viable to way to use computational re- 
sources more efficiently. This result was verified through 
simulation of a Schwarzschild black hole in geodesic co- 
ordinates, for which we have an analytic solution for 
comparision, and through simulations of Schwarzschild 
in a variation of the 1-blog (singularity avoiding) slic- 
ing. By nesting several levels of mesh refinement regions, 
we are able to convergently resolve the puncture while 
simultaneously pushing the boundary of our domain to 
120 M and keeping the computational problem size mod- 
est. Combined with our earlier results showing that grav- 
itational waves pass through such FMR interfaces with- 
out significant reflections [? ], we have now studied, in 
detail, the effects of FMR interfaces on the two primary 
features, waves and strong potentials, of astrophysically 
interesting spacetimes. 

Our method for handing the interface condition, based 
in part on the Paramesh infrastructure, is detailed. We 
find that, in handling the interface condition between 


FMR levels, it is sufficient to interpolate with one order 
lower accuracy than would naively be expected without 
spoiling the convergence of our code. This leads to the 
somewhat counter-intuitive result that derivatives of our 
evolved fields appear non-convergent at grid points ad- 
jacent to the interface (c.f. Figs. 3, 5, 1-blog), but that 
the evolved variables, for which these same derivatives 
are sources, remain convergent at all points in the inte- 
rior of our domain (c.f. Figs 2, 5, and 1-flog). Because it 
is the fields themselves, not their derivatives, that repre- 
sent the physical solution, we find no problems with this 
understood behavior. 

In the course of our investigation, we have also ob- 
served what we call the “point- two M problem” in 1-flog 
runs. At intermediate times (GIVE AN ESTIMATE), 
high frequency noise appears at radial coordinate r w 
0.2 M in the lapse a and the trace of the extrinsic curva- 
ture K. This noise then dissipates by later times so that, 
if the data is viewed outside of a certain window of time, 
it will not be observed. Higher resolution exacerbates 
this problem. We have verified that this phenomenon 
is present both with and without mesh refinement in- 
terfaces in the domain, and have further verified that 
it appears in an independently written, one-dimensional 
code. We therefore conclude that this is a real feature 
of the gauge choice used. Having chosen a generally ac- 
cepted gauge, and having focused on effects of the mesh 
refinement interfaces in this work, we have not fully inves- 
tigated the cause of nor possible remedies for this appar- 
ent pathology. We note it here, however, as an interesting 
topic for future investigation. 
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